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Основываясь на методах решения некорректных задач, восстановлены теплофизические и массообменные характеристики про- 
мерзающих-протаивающих дисперсных грунтов, разработаны алгоритмы, которые использованы для решения прямых задач 
тепломассообмена. 


Введение 

Освоение северных территорий ведется с техно- 
генным загрязнением (промстоками, нефтепро- 
дуктами, радионуклидами аварийного подземного 
ядерного взрыва и другими экологически опасны- 
ми загрязнителями) мерзлых грунтов. При гло- 
бальном потеплении климата указанные процессы 
усиливают деградацию многолетней мерзлоты, 
приводят к заболачиванию огромных территорий, 
потере несущей способности грунта, миграции 
экологически опасных загрязнителей в речную 
систему и повышению вероятности попадания их в 
организм человека. Принятие различных инженер- 
но-технических решений, меры по борьбе с загряз- 
нением требуют усовершенствования методологии 
численного математического моделирования с уче- 
том реального физического процесса промерзания 
- протаивания порового раствора грунта. 


1. Постановка задачи 


В механике мерзлых грунтов широкое практичес- 
кое распространение имеют математические модели 
многофазных сред [1]. Для получения математичес- 
кой модели тепломассообмена применяют законы 
сохранения массы, импульса и энергии для каждой 
фазы, а затем, суммируя, получают общее уравнение 
сохранения для всей системы. При этом систему 
уравнений выписывают в рамках следующих предпо- 
ложений [2]: скелет пористой среды - лед, поровая 
жидкость и воздух (газ) слабосжимаемы; пренебрега- 
ют связанной водой; при образовании льда примесь 
отторгается в объеме незамерзшей воды и в осадок не 
выпадает. Система уравнений замыкается равновес- 
ным термодинамическим условием существования 
льда, порового раствора и воздуха. В мерзлых диспе- 
рсных грунтах большую роль играет физическое сос- 
тояние связанной воды, которое усиливается в зави- 
симости от типа техногенного загрязнения. В зоне аэ- 
рации (в деятельном слое) обычно пренебрегают воз- 
духом и сжимаемостью пористой среды (скелета), ль- 
да и воды. Тогда из модели многофазной среды выте- 
кает общеизвестная модель А.В. Лыкова. 
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О < т < Н, т т > т > 0. 

Уравнения (1-3) замыкаются равновесным 
уравнением незамерзшей воды: 

К=К в =К е (тл,С) (4) 


и условием замерзания соли вместе с почвенной влагой: 

С,=к. ІХ С г (5) 

где С и С е - объемная теплоемкость грунта и воды, 
Дж/(м 3 -К); Т - температура, К; Я - теплопровод- 
ность, Вт/(м-К); г - время, с; X - теплота фазового 
перехода, Дж/м 3 ; ]Ѵ - суммарная влажность, % 
(масса всей влаги/ масса скелета); Ѳ - объемная 
влажность, % (отношение объема влаги к объему 
пор); Н=Р-і- напор, м; Р - всасывающее давление 
влаги, м; г, 7, (ось г направлена вертикально вниз) - 
пространственные координаты, м; к - коэффици- 
ент диффузии, м 2 /с; к ф - коэффициент фильтрации, 
м/с; Р - коэффициент конвективной диффузии 
примеси, м 2 /с; Ѵ=(Ѵ Г ,К) - скорость фильтрации, 
м/с; С„ С, - концентрации примеси в воде и льду, %; 
N - концентрация примеси в твердой фазе, %; 
Р - коэффициент скорости обмена, 1/с; к л (ХД - 
(межфазный) коэффициент распределения примеси. 

Уравнение (1) учитывает процесс промерзания - 
протаивания порового раствора с учетом фильтра- 
ции жидкой фазы. Движение самого порового раст- 
вора (воды) с учетом льдовыделения описывается 
аналогичным уравнением параболического типа 
(2). Выражения (3), (5) характеризуют процесс со- 
лепереноса в промерзающих - протаивающих грун- 
тах. Выражение (2) называется потенциальной, (2*) 

- влажностной формами уравнения влагопереноса. 
Потенциальная форма (2) получает из влажностной 
(2*) при замене переменной. Для потенциальной 
формы нужно установить функциональные зависи- 
мости объемной влажности и коэффициента 
фильтрации от давления, а для влажностной формы 

- зависимость от суммарной влажности, которая 
легко восстанавливается на практике. 
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Ключевыми моментами при построении мате- 
матических моделей тепломассопереноса в про- 
мерзающих и протаивающих грунтах являются: 
способ локализации области фазового перехода; 
выбор формы представления уравнения влагопере- 
носа. Известны два ( структурных) подхода: соглас- 
но первому из них фазовый переход локализован 
на поверхности раздела фаз (при определенной 
температуре); второму - фазовое превращение 
происходит в протяженной области (так называе- 
мая модель фазового перехода в спектре темпера- 
тур). Становление второго подхода связано в пер- 
вую очередь с экспериментальными работами, ко- 
торые свидетельствовали о протекании этих про- 
цессов в мерзлых грунтах. А. Г. Колесниковым [3] 
предложена математическая формулировка задачи 
о температурном режиме грунта, промерзающего в 
спектре температур. В дальнейшем температурная 
задача была дополнена уравнением массопереноса. 

Вначале рассмотрим, какие трудности встреча- 
ются при различных видах (структурных подходах) 
задания математических моделей теплосолевлаго- 
переноса. При потенциальной форме записи урав- 
нения влагопереноса возникает проблема представ- 
ления его в зоне промерзания - протаивания. 
Влажностная форма математической модели учи- 
тывает процесс влагосолепереноса только в талой 
зоне и поэтому трудно оценить, сколько воды пере- 
ходит в лед (или, наоборот) при промерзании (про- 
таивании). А также неизвестно, сколько соли захва- 
тывается льдом при промерзании порового раство- 
ра грунта. В связи с этим, возникает "проблема за- 
дания граничного условия" для влаги и соли на под- 
вижном фронте фазового перехода со стороны та- 
лой зоны [4]. Для решения задач теплосолемассооб- 
мена в насыщенных средах В.М. Ентов, А.М. Мак- 
симов, Г. Г. Цыпкин [5] успешно применили второй 
подход (с зоной промерзания - протаивания). При 
этом проблема граничных условий первого подхода 
автоматически снимается, но вместе с тем возника- 
ет проблема неопределенности массообменных ха- 
рактеристик в зоне промерзания - протаивания. В 
связи с этим они с большим удивлением писали 
"почему все это не было детально изучено ранее?”. 
Причиной этого считаем отсутствие соответствую- 
щих приборов и установок для экспериментального 
изучения и методов идентификации параметров 
модели. Адекватность вышеуказанных моделей 
проверена в лабораторных и натурных условиях и 
имеет хорошую сопоставимость [6]. Параметры вы- 
ражений (4) и (5) находятся из эксперимента. 

2. Идентификация параметров модели 

Восстановление теплофизических и массообмен- 
ных характеристик с учетом процесса промерзания - 
протаивания порового раствора относится к классу 
нелинейных задач с сильноменяющимися коэффи- 
циентами. Традиционными методами их восстанов- 
ления невозможно, поэтому воспользовались 
экстремальными методами минимизации функцио- 
нала [7]. При фазовом переходе порового раствора 


все теплофизические (теплоемкость, теплопровод- 
ность) и массообменные (коэффициенты диффузии, 
фильтрации и конвективной диффузии примеси) ха- 
рактеристики выражаются через функцию количест- 
ва незамерзшей воды. Для простоты изложения рас- 
смотрим процесс идентификации параметров функ- 
ции количества незамерзшей воды от температуры. 


Как известно, процесс протаивания (промерза- 
ния) влажной мерзлой (талой) дисперсной среды 
сопровождается переносом влаги, которая приво- 
дит к существенным изменениям теплофизических 
характеристик. При проведении эксперимента пе- 
реносом влаги можно пренебречь, если будут соб- 
людены следующие условия: 

- Во-первых, эксперимент проводится в цилиндри- 
ческой ячейке (сосуде) достаточно малого радиуса 
(Д<0,015 м, высота к> 6К.), такого, что поток вла- 
ги, мигрирующий по радиусу, незначителен. 

- Во-вторых, создается равномерный, достаточно 
интенсивный тепловой поток на поверхности из- 
мерительной ячейки. В этом случае влага не успева- 
ет мигрировать за время проведения эксперимента 
ввиду его сравнительной кратковременности. 

- В-третьих, при проведении эксперимента 
должны быть соблюдены условия теоремы 
единственности решения обратной задачи теп- 
лопроводности Н.В. Музылева [8]. 
Температурное поле исследуемого образца в ци- 
линдрическом сосуде описывается следующим 
уравнением теплопроводности: 



0 < т < т т , 0 <г < К, (6) 

0< г <т т , (7) 


-г ѵ А(г)^ Г =г ѵ ^г(г) 
ог 

Г=К 


0 < г < г 


( 8 ) 


Г(г,0)= Г°(г) 0 <г<К, (9) 


где 

сСП = Сэф =с ск +с,Г +( С „ -с г )Ж ж (Г)+1^-, (Ю) 


Л(Т) = Л м +(Л т -Л м ) 


и’°аи; г 


(П) 


ѵ = 0, 1 - соответствует пластине и цилиндру. 

На границе при г = 0 выполняется условие огра- 
ниченности решения (7). 

Требуется восстановить параметры 
(с м ,с т ХЛт,Ке(Т))Т=(Щ,Щ,---,и щ У =и по известным 
замерам температуры Т°=Т}(г,т) в точках г і (/=1 ,п т ). 

Задачу (6-9) сформулируем в виде задачи опти- 
мального управления. В качестве оптимальности 
выберем целевой функционал среднеквадратичес- 
кого отклонения замеренных температур от расчет- 
ных значений: 
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^ (и) = ^ I р, (г (и)- Т; ) Ат = Ш- Т э || , «е[/с й”°, (12) 

где р, - весовые множители; 7Хи), 7] э - расчетные и 
замеренные температуры; число замеров по ра- 
диусу образца; Се і" г [0,г т ] - пространство замеров. 

Разработана следующая последовательность отыс- 
кания неизвестных: из термограммы опыта выбира- 
ются области мерзлых и талых зон и в каждой из них 
определяются коэффициенты (с„Д„) и (с г Д г ); затем с 
помощью (10) вычисляются с ск ,ѴС пс \ на последнем эта- 
пе, используя все замеры температуры, восстанавли- 
ваем параметры функции незамерзшей воды. 

Минимизацию целевого функционала можно 
осуществить методами скорейшего спуска, сопря- 
женных градиентов или модифицированным мето- 
дом скорейшего спуска [7,9]. 

График восстановленной функции количества 
незамерзшей воды лежит между кривыми 2, 3 
(рис. 1), полученными в результате расчета количест- 
ва незамерзшей воды, начиная с талой (баланс воды) 
и мерзлой (баланс льда) зон, соответственно. Неко- 
торые различия между ними наблюдаются в области 
начала и конца фазовых переходов. Расхождение в 
начале фазовых переходов можно объяснить следую- 
щим образом. При температуре 273 К замерзает зна- 
чительная часть свободной воды и кривая незамерз- 
шей воды имеет "скачкообразный" вид. Этот "ска- 
чок" при дискретной реализации приближенно ап- 
проксимируется дельта функцией (линейно). 



Рис. 1. Сравнение результатов восстановления функции ко- 
личества незамерзшей воды, полученным числен- 
ным (1) и квазистационарным (2,3) методами 


ными, но меняется вид функции плотности расп- 
ределения источника тепла. 

Мощность источника тепла выражается условием: 

|<5(Г, ду/Г = 1 или 1 1 с ІѴнв( ' Т,С) АТ = 1. 


В реальных условиях плотность распределения 
зависит от дисперсности, минерального состава 
грунтов, типа растворенных ионов и т.д. В резуль- 
тате обработки данных экспериментальных иссле- 
дований при засолении песка хлоридом натрия по- 
лучены различные виды графиков плотности расп- 
ределения источника тепла (рис. 2). С повышением 
концентрации растворенных солей усиливаются 
растворение льда и понижение средней температу- 
ры фазового перехода. Идентификация параметров 
функции количества незамерзшей воды не занима- 
ет много времени и труда [9] . 


Аппроксимационные формулы (5-функции со- 
ответственно имеют вид (рис. 2): 
а) 

[— , |г|<д, 

8(Т,А)= 2Д 1 1 

[о, |г| > д, 
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т 4 <т, 


В вычислительной теплопередаче широкое 
практическое применение получил метод сглажи- 
вания [10, 11]. Это равносильно предположению, 
что фазовый переход происходит не при одной оп- 
ределенной температуре Т ф , а в некотором интерва- 
ле температур, длина которого определяется вели- 
чиной параметра сглаживания А. При численном 
эксперименте обычно задают различные значения 
этого параметра, при этом изменяется ширина ди- 
апазона температур фазового перехода. С другой 
стороны, мощность и средняя температура источ- 
ника тепла фазового перехода остаются постоян- 


где А, а, Ъ, А, А,, 7] (/=1,2, 3,4), Щ,, Щ 1С - параметры 
функции количества незамерзшей воды. 

Из рис. 2 видно, что с повышением концентра- 
ции однородного раствора (хлорида натрия) содер- 
жание незамерзшей воды увеличивается и её за- 
мерзание происходит при температуре эвтектики. 
Идентифицированные расчетные формулы учиты- 
вают реальный процесс промерзания - протаива- 
ния порового раствора мерзлого грунта в спектре 
температур. Поэтому предложенный способ можно 
назвать естественным методом сглаживания. 
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Рис. 2. Плотность распределения источника тепла: а) прямоугольник; Ь) парабола; с) С, = 5; 6) 15; е) 20; {) 30 % 


Сама граница фазового перехода явно не выде- 
ляется и не участвует в построении разностной схе- 
мы, а при необходимости свободная граница иден- 
тифицируется как изотерма средней температуры 
фазового перехода Тф. 

с 1 1 ^ ё\Ѵ нв (Т, С) ігті 

Т с , = г т — НвУ (1 Т 

Ф Ш -Ш 1 яг 

"о ГГ ПС 0 и1 

3. Численная реализация 

Для численного решения системы дифференци- 
альных уравнений (1-3) в области О=[0,Н]х[0,т т ] 
введем неравномерную сетку: 


СО/ 1т = СО /, X С0 Т , Ю/, = {х 0 = 0, х,- = Х;_\ +1ц, І=\,Щ, 
со т = {г у =ут, у =1,2,3,...}, %і ={1ц +* - )/2, «■= 1, N 

На множестве внутренних узлов щ={х\, і= 1 ,7Ѵ— 1} 
систему уравнений аппроксимируем разностной 
схемой [12]: 

С эф Т т +В 1 Т +с в С 1 Т = (р 1 , 

5+1 5 5 V 

ц^-^- + ^-^ + І) 2 Р + {кф), =<р 2 , 

(іу т+ о;іу= ( р 2 ) 

8 Т + ГХБ + С 2 8 = <р, 
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где 

ОуТ =-(АГ-)_ ;> ДЖ = -(к фѢ (Р),)., 

АЖ = -(АА,Ж)А, А 5 = -ФААХ-);, 
с,г = (УХ+Ѵ+Т-), СА = ((ЖТ7А)_- + АААХА 
ср { = ищ, <р 2 = % = о, 

5 = (ж в с в +ж л с л )/іоо, ж = ж і+ ж, ѳ = ѳ +ѳ, 

К=КЛП С в =і і2 5,Ѵ = Ѵ + + Ѵ-, 

Ѵ + =0,5 (Ѵ+\Ѵ\)>0, Ѵ~ =0,5(Ѵ-\Ѵ\)<0, 

і = % , 7?, = 1 - і, ъ = 1 00 /(АЖ, + Ж в ). 

Разностную схему температурной задачи получили 
следующим образом. Используя функцию незамерз- 
шей воды Ж в =Ж^,( Г, Ж,5) и учитывая независимость 
потенциалов Г, Ж 5, множитель теплоты фазового пе- 
рехода д\Ѵ,/дт ур. (1) можно записать в виде: 

дЖ, _ ЗЖ 7 д Т ЗЖ, діу д!У л 

дт дТ г ЗЖ дт 35 дт ' 

Отсюда при помощи частных производных 
функции Ж/ 

ЗЖ л _5(Ж-Ж в )_ ЗЖ в 
~дТ~ дТ “ ~дТ' 

зж,_ з»(Г)ж 

ЗЖ ЗЖ 1 Л 

зж л _5(ж-ж в )_ зж в 
“з/ _ дТ^~~~дЗ~ 

получим: 

5 ж, зжзг . зж зж^5 

Зг ЗГ Зг ѵ ■ Зг 35 Зг ' 

Первое слагаемое данного выражения учитыва- 
ется в эффективной теплоемкости (10). Последним 
членом выражения можно пренебречь, т.к. его вклад 
в реальных условиях достаточно мал (д]Ѵ нв /д5~0). 

В левой части уравнения фильтрации (2) при- 
меняется метод Ньютона для локализации кривой 
водоудерживания (/и=дѲ/дР). Введение 3-функции 
для приближенного решения задачи Стефана поз- 
воляет строить разностные схемы со сглаженными 
коэффициентами, т.е. совершается переход к 
обычной задаче параболического типа. Преимуще- 
ством данного подхода является то, что не нужно 
выбирать параметр сглаживания А. Численная реа- 
лизация нелинейной задачи (1-3) осуществляется с 
помощью итерационных методов сквозного счета, 
что намного упрощает процесс решения многомер- 
ных задач тепломассообмена в мерзлых грунтах 
при использовании экономичных аддитивных ло- 
кально-одномерных разностных схем. Через функ- 
цию количества незамерзшей воды выражаются 
доли порового раствора ц,( Т) и засоленности гі 2 ( 7), 
которые участвуют в процессе фильтрации (диф- 
фузии). Разностные операторы уравнений диффу- 


зии и солепереноса имеют диагональное преобла- 
дание по столбцам, и достаточные условия устой- 
чивости выполняются в гильбертовом простран- 
стве АА), А и/) [13]. 

Динамика температуры 

0,0 5,0 10,0 г 15,0 20,0 25,0 




Рис. 3. Распределения температуры и суммарной влажности 
при С в =0 

Динамика температуры 


0,0 5,0 10,0 Г 15,0 20,0 25,0 



О 0,0 5,0 10,0 г 15,0 20,0 25,0 



0,0 5,0 10,0 г 15,0 20,0 25,0 



Рис. 4. Распределения температуры, суммарной влажности 
и засоленности при С в =100 г/л 
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Для проверки возможностей математической 
модели приведем пример. Основная цель - мето- 
дическая: показать какого рода численные иссле- 
дования подземных фильтрационных процессов 
можно проводить в промерзающих - протаиваю- 
щих грунтах и насколько эффективно работает при 
этом предложенный алгоритм. Алгоритм решения 
при промерзании-протаивании мерзлого грунта 
построен с учетом знака скорости фильтрации с 
помощью направленных разностей. Численные 
расчеты проведены с учетом природно-климати- 
ческих условий, радиационно-теплового баланса 
поверхности и массообменных характеристик 
грунтов Центральной Якутии [6]. 

Рассмотрим случай, когда в грунт поступает (в 
летний период) надмерзлотная загрязненная грун- 
товая вода при различных концентрациях. Техно- 
генная вода, загрязненная промстоками, проника- 
ет по полосе 0, 5. ..0, 75 м с напором - 0,5 м и фильт- 


СПИСОК ЛИТЕРАТУРЫ 

1. Нигматулин Р.И. Основы механики гетерогенных сред. — М.: 
Наука, 1978. —336 с. 

2. Васильев В.И., Максимов А.М., Петров Е.Е., Цыпкин Г.Г. Теп- 
ломассоперенос в промерзающих и протаивающих грунтах. — 
М.: Наука, 1997. -224 с. 

3. Колесников А.Г. К изменению математической формулировки 
задачи о промерзании грунта // Доклады АН СССР. — 1952. — 
Т. 32. -№6. -С. 889-891. 

4. Гречищев С.Е., Чистотинов Л. В., Шур Ю.Л. Криогенные фи- 
зико-геологические процессы и их прогноз. — М.: Недра, 1980. 
-284 с. 

5. Ентов В.М., Максимов А.М., Цыпкин Г.Г. Образование двух- 
фазной зоны при промерзании пористой среды. — М.: ИПМ 
АН СССР. Препринт № 269, 1986. -56 с. 

6. Пермяков П.П., Романов П.Г. Тепло- и солеперенос в мерзлых 
ненасыщенных грунтах. — Якутск: ЯФ Изд-ва СО РАН, 2000. 
— 126 с. 

7. Алифанов О.М., Артюхин Е.А., Румянцев С.В. Экстремальные 
методы решения некорректных задач. — М.: Наука, 1988. —288 с. 


рирует вниз в сторону реки. На рис. 4 приведены 
распределения температуры, суммарной влажнос- 
ти и засоленности грунта (кг/м 3 ) в конце июля при 
коэффициенте фильтрации к ф = 0,48 м/сут. С повы- 
шением концентрации промстока процесс раство- 
рения льда усиливается, что сопровождается пони- 
жением температуры, а также увеличением под- 
вижности жидкой фазы. Численный эксперимент 
проведен с учетом сезонной динамики уровня ре- 
ки. Из рисунка видно, что поступление техноген- 
но-загрязненной воды действует на динамику тем- 
пературного и влажностного режимов массива. 

Вывод 

Предложенные алгоритмы решения и програм- 
мные средства можно применять для прогноза дина- 
мики техногенного загрязнения в мерзлых грунтах. 
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